Phylogenomic insight into dysploidy, speciation, and plastome evolution of a small Mediterranean genus Reichardia (Cichorieae; Asteraceae)

Reichardia Roth is a small Mediterranean genus comprising ten homogeneous species with basic chromosome numbers of 7, 8, and 9. To assess the plastid genome evolution and differentiation of Reichardia species, we assembled the complete plastome sequences of seven Reichardia and two Launaea species and conducted various phylogenomic analyses comparatively with nuclear ribosomal DNA ITS sequences. Reichardia and Launaea plastomes were highly conserved in gene content and order, containing 130 genes. Plastid phylogenomic reconstruction strongly suggested that Reichardia was a sister to Launaea, and its common ancestor initially diverged into two major lineages: the first containing species with n = 8 chromosomes exclusively, and the other with n = 9, 8, and 7 chromosomes. Although the ancestral Reichardia karyotype was suggested to most likely be n = 9 from ancestral chromosome number reconstruction, the pattern of descending dysploidy indicated by the phylogenetic trees based on nuclear ribosomal DNA ITS was less evident in the trees based on the plastome. Possible reasons for these findings are discussed.

www.nature.com/scientificreports/ been investigated within a molecular phylogenetic framework. Siljak-Yakovlev et al. 19 plotted the karyological, cytogenetic, and pollen traits of Reichardia species onto the nuclear ribosomal DNA (nrDNA) internal transcribed spacer (ITS) phylogeny and suggested that the descending dysploidy in Reichardia chromosomal numbers was accompanied by changes in heterochromatin patterns and modifications of the location and organization of ribosomal genes. They also reported smaller pollen sizes in R. picroides and R. intermedia, which correlated with their smaller genome sizes and the reduced chromosome number within the genus (n = 7 for both species); however, in our view, these reports for R. intermedia may be erroneous. Moreover, these cytological studies did not include appropriate closest relatives or progenitors of the genus. Without robust phylogenetic relationships between Reichardia and other closely related genera and precise species relationships within the genus, it may be unsuitable to explore genome evolution and differentiation. The phylogenetic position of Reichardia was previously inferred using nrDNA ITS and chloroplast DNA (cpDNA) psbA-trnH intergenic spacer and matK gene sequences [25][26][27][28] . In all phylogenies, the genera Launaea and Reichardia represented early diverging lineages within the subtribe Hyoseridinae, and Reichardia was monophyletic, which is concordant with the first comprehensive molecular analyses of the tribe Cichorieae 29 . The ITS tree strongly supported the monophyly of Reichardia 25,26,28 , which was sister to the clade containing the genus Launaea and the remaining genera of the subtribe Hyoseridinae, that is, the currently circumscribed Sonchus sensu lato (s.l.). The cpDNA trees based on either the psbA-trnH intergenic spacer 27 or matK 28 gene regions also supported the monophyly of Reichardia and the basal position of Reichardia and Launaea. However, cpDNA phylogenies in previous studies 27,28 were too limited to determine the phylogenetic position of Reichardia and robust species relationships due to insufficient molecular markers and species representation. The relationships between species and genera were either highly unresolved or weakly supported, specifically in the psbA-trnH tree 27 . Unlike the nrDNA ITS-based phylogeny, Reichardia and Launaea shared the most recent common ancestor, and these two genera were sister to the remaining genera of the subtribe Hyoseridinae, which was also weakly supported (< 50% in the psbA-trnH tree) 27 or moderately supported (86% in the matK tree) 28 . These conclusions are consistent with subsequent studies based on ITS datasets 19,30 .
With the advent of high-throughput sequencing technologies for next-generation sequencing (NGS), massive amounts of data have now become available based on the considerable genome-wide variation of entire plastid genomes through plastome sequencing. The benefits of genome-wide variation could increase phylogenetic resolution tremendously and significantly enhance our understanding of plant evolution and diversity in the field of plastid genetics and genomics 31 . Comparative genomic analysis using whole plastome sequences is now an efficient option to improve phylogenetic resolution at lower taxonomic levels that are currently hindered by limited sequence variation due to recent divergence, rapid radiation, and conservative genome evolution of plastomes 32 .
In this study, we sequenced and assembled the whole plastid genomes of seven Reichardia species in addition to R. ligulata, sequenced in our previous study 33 . The plastomes of two additional species from the closely related genus Launaea were also sequenced [25][26][27][28] . Based on these complete plastome sequences, we first tested the previous phylogenetic hypotheses proposed by nrDNA ITS and cpDNA markers; plastome-based phylogeny was reconstructed using genome-wide variation and compared with the nrDNA ITS phylogeny, newly reconstructed in this study. Second, through the ancestral karyotype reconstruction, we inferred the karyotype evolution in the dysploidy process of Reichardia species in the framework of the plastome and nrDNA ITS-based phylogenies. We also performed several comparative plastome analyses to determine the structure, gene content, and rearrangements of the plastid genomes of Reichardia and two Launaea species. In addition, we identified highly variable plastid regions that could be utilized as useful markers for further population genetic or phylogeographic studies of Reichardia and closely related genera.

Results
Genome features, content, order, and organization of Reichardia plastomes. The eight plastomes of Reichardia species (R. albanica, R. crystallina, R. famarae, R. gaditana, R. intermedia, R. ligulata, R. picroides, and R. tingitana) and two species of Launaea (L. arborescens and L. nudicaulis) were highly conserved in terms of gene content and arrangement, with 98.6% pairwise sequence similarity (99.2% among eight Reichardia species only) ( Table 1). The total length of eight Reichardia plastomes ranged from 152,067 (R. albanica) to 152,652 (R. famarae) base pairs (bp) with an overall guanine-cytosine (GC) content of 37.6%, whereas two plastomes of Launaea were slightly longer than Reichardia plastomes, that is, 153,119 bp for L. nudicaulis and 154,086 bp for L. arborescens with an overall GC content of 37.4%. The plastomes of Reichardia and Launaea species consisted of four typical regions: LSC, SSC, and a pair of inverted repeat regions (IRs), sharing exactly the same genes and similar gene contents at all boundaries among the four regions, with slight changes in the length of intergenic regions. They all contained the functional protein-coding gene ycf1 at SSC/IR with its pseudogene copy, ycf1Ψ at IR/SSC, and functional rps19 at LSC/IR with pseudogene copy rps19Ψ at IR/LSC endpoints (Fig. 2).
Each of the eight Reichardia and two Launaea cp genomes contained 130 genes, including 87 protein-coding genes (excluding pseudogenes), six rRNA genes, and 37 tRNA genes (Tables 1 and 2). Twenty-four genes contained introns, including nine tRNA genes. Three genes, clpP, rps12, and ycf3, had two introns. The trnK-UUU tRNA gene harbors the largest intron, which contains the matK gene. In total, 17 genes were duplicated in the IR regions, including seven tRNAs, three rRNAs, and seven protein coding genes. The trans-splicing gene rps12, consisting of three exons, was located in the LSC region of exon 1, whereas exons 2 and 3 of the gene were embedded in the IR regions. Parts of ycf1 and rps19 duplicated in the IR region were considered pseudogenes in all cp genomes sequenced in this study.  Table S1). The codon usage numbers in the ten Reichardia and Launaea plastomes ranged from 22,765 to 22,791, and the patterns of frequently used codons were also consistent among them except for slight variations in codon usage numbers. Within Reichardia, species with chromosome numbers of n = 8, except R. intermedia, showed lower codon usage (22,765) than species with chromosome numbers of n = 9 (R. albanica: 22,781) or n = 7 (R. picroides: 22,790), and R. intermedia (22,791). Codon usage values are described by relative synonymous codon usage (RSCU), which reflects how often a particular codon is used relative to the expected number of times that codon would be used in the absence of codon usage bias. All the RSCU values for each amino acid considered in this study were similar among the ten plastomes. The highest RSCU value was indicated in the usage of the UUA codon for leucine (1.9-1.92), followed by AGA for arginine (1.84-1.86), while the lowest were indicated in the usages of the AGC codon for serine (0.33-0.36) and CUG for leucine (0.35-0.38). The possible RNA editing sites predicted among the Reichardia and Launaea plastomes ranged from 43 to 48 in the 19 genes of the 35 protein-coding genes ( Fig. 4 and Supplementary Table S2). We found that the RNA editing patterns across the ten Reichardia and Launaea plastomes were similar in gene location and codon conversion type of the predicted RNA editing sites, and only slight changes were observed in the number of editing sites for several codon conversions. Within Reichardia, all species with chromosome numbers of n = 8 showed lower numbers of predicted editing sites (43 for all species, with the exception of 44 for R. intermedia) than species with chromosome numbers of n = 7 (45 for R. picroides) and n = 9 (R. albanica), with 47 sites showing the highest number in the genus. Genes with predicted sites included photosynthesis-related genes (atpA, ndhA, ndhB, ndhD, ndhF, ndhG, petB, psbF, and psbL), self-replication genes (rpl20, rpoA, rpoB, rpoC1, rpoC2, rps2, and rps14), and others (accD, ccsA, and matK). We detected no RNA editing sites in atpB, atpF, atpI, clpP, petD, petG, petL, psaB, psaI, psbB, psbE, rpl2, rpl23, rps8, rps16, and ycf3. The highest numbers of potential editing sites were found in the NADH dehydrogenase genes, which was consistent with previous findings in tobacco, maize, rice, and other plants [34][35][36][37] ; the ndhB gene was the highest at 9-10 sites, followed by the ndhD gene at 4-5 sites. The highest conversions in the editing frequencies of codons associated with the corresponding amino acid changes were represented by the changes from serine (S) to leucine (L) (average confidence score of 14.52), followed by proline (P) to leucine (L) (average confidence score of 8.76).
The divergence level of nucleotide diversity between the Reichardia and Launaea plastomes was compared using DnaSP 38 and visualized by plotting with mVISTA 39 . The results showed a high degree of synteny and gene order conservation in the mVISTA graph (Fig. 5). The overall nucleotide diversity (Pi) of the ten plastomes was 0.00520, with 2586 polymorphic sites, ranging from 0 to 0.02836, which was higher than that of Sonchus species Table 1. Summary of the genomic characteristics of two Launaea and eight Reichardia chloroplast genomes used for comparative genomic analyses in this study. Herbarium code: APP: Parco Nazionale del Gran Sacco e Mont della Laga-Universita di Camerino, Italy, OS: Ohio State University herbarium, SEV: Herbario de la Universidad de Sevilla, TEX: University of Texas Herbarium. a Sequenced from a previous study 33 . The remaining nine cp genomes were sequenced for this study.  www.nature.com/scientificreports/   www.nature.com/scientificreports/ (average Pi value of 0.00283, ranging from 0 to 0.01593) belonging to the same subtribe Hyoseridinae 33 . Genetic polymorphisms in different regions of the chloroplast genome varied substantially. The SSC region, where the most variable gene, ycf1, was located, was the most divergent (average Pi, 0.0102), whereas the two IR regions were highly conserved (average Pi, 0.00173). Six divergence hotspots among the Reichardia and Launaea plastomes are suggested as potential chloroplast markers: five intergenic regions (trnH-psbA, trnC-petN, ndhC-trnV, trnL-rpl32, and rpl32-ndhF) and one protein-coding gene (ycf1) (Fig. 6).  www.nature.com/scientificreports/ Selective pressure in genes or genomic regions is inferred by the proportion of amino acid substitutions, and the ratio (ω = dN/dS) of nonsynonymous substitution (dN) and synonymous substitution rates (dS) has been widely used as a genomic signature of selective pressure acting on a protein-coding gene, that is, ω = 1 indicates neutral mutations, ω < 1 indicates purifying selection, and ω > 1 indicates diversifying positive selection 40 . We identified that seven genes potentially evolved under positive selection in ten Reichardia and Launaea plastomes by calculating the dN/dS ratio using various site-specific substitution models implemented in EasyCodeML (Table 3 and Supplementary Table S3) 41,42 . These genes included NADH-dehydrogenase subunit genes (ndhB and ndhG), subunit genes of photosystem I (psaI) and II (psbH), and large subunit genes of ribosomal protein (rpl20), ycf2, and ycf15. Positively selected sites were suggested based on the posterior probability calculated using the Bayes empirical Bayes (BEB) method 43 with a cutoff > 0.95 and > 0.99 indicated with asterisks (* and **, respectively) in Table 3. Despite the critical importance of the genes that a plastome carries and its high conservativeness, the latest empirical evidence revealed that the variable genes potentially evolve under positive selection in the plastomes of a few other plant groups; three genes (rps2, rbcL, and ndhG) have been identified in Paulownia 44 , five (rbcL, clpP, atpF, ycf1, and ycf2) in Panax 45 , three (clpP, ycf1, and ycf2) in the tribe Sileneae 46 , and six (accD, rbcL, rps3, ndhB, ndhD, and ndhF) in Rosaceae 37 .
Phylogenomic analyses and ancestral state reconstructions. Phylogenetic trees were constructed using maximum likelihood (ML) and Bayesian inference (BI) methods based on the plastome and nrDNA ITS sequences to infer the phylogenetic relationships among Reichardia and the closely related species in the subtribe Hyoseridinae, with Lactuca sativa as an outgroup. The topological structures of ML phylogenetic trees constructed by IQ-TREE 47 and BI trees by MrBayes 3.2.7a 48 were consistent, and BI trees showed higher support values than ML trees ( Fig. 7 and Supplementary Fig. S1). The best-fit evolutionary models in ML analyses were www.nature.com/scientificreports/ selected as "TVM + F + I" for plastome trees and "TIM3e + G4" for the nrDNA ITS tree using ModelFinder 49 implemented in IQ-TREE. The plastome phylogenies were reconstructed using both whole plastome sequences (Fig. 7a) and 80 chloroplast protein-coding genes ( Supplementary Fig. S1) of 18 representative Hyoseridinae species. The phylogenies based on plastome datasets shared the same topology and provided much greater resolution with strong bootstrap (BS) values for inter-and intra-generic relationships compared to previous studies 27, 28 . The genus Reichardia was resolved as monophyletic and shared the most recent common ancestor with the genus Launaea, placing the clade of Reichardia and Launaea sister to the genus Sonchus s.l. Within Reichardia, R. albanica (n = 9) shared the most recent common ancestor with R. picroides and R. intermedia with chromosomal numbers of n = 7 and 8 robustly (99% or 97% BS support value in ML trees and posterior probability, PP, 1 in BI trees). The remaining five species (R. ligulata, R. famarae, R. crystallina, R. tingitana, and R. gaditana) with chromosome numbers of n = 8 formed a clade with full support (100% BS for ML and PP 1 for BI phylogenies), which was sister to the clade of n = 7, 8, and 9. On the other hand, the nrDNA ITS phylogeny showed incongruence in inter-and intra-generic relationships compared with the plastome phylogeny (Fig. 7b). Unlike plastome phylogenies, Launaea shared the most recent common ancestor with Sonchus s.l., instead of Reichardia; this relationship, however, was weakly supported (68% BS in ML and 0.7867 PP in BI trees). Reichardia was monophyletic in both nrDNA ITS ML and BI phylogenies (100% BS and 1 PP, respectively), but differed  The ancestral chromosome numbers optimized on the plastome and nrDNA ITS ML phylogenies also reflected the differences between both phylogenies. The ancestral chromosome number of the most recent common ancestor (MRCA) shared by three genera, Reichardia, Launaea, and Sonchus, was inferred to be equivocal with 2n = 14, 16, 18, and 32, and 2n = 36 for both plastome (node P21) and nrDNA ITS (node 29) datasets (Fig. 8). Within Reichardia, ancestral state reconstruction using the nrDNA ITS tree suggested that the MRCA of Reichardia most likely had 2n = 18 (56.5% probability at node 40), suggesting a decrease in dysploidy (Fig. 8b

Discussion
This study provides, for the first time, greater resolution with high support for the plastome-based inter-and intra-generic relationships among Reichardia and closely related genera. All three genera, Reichardia, Launaea, and Sonchus, were confirmed as monophyletic with strong support values (BS 100% each) in both phylogenies of full plastome sequences and concatenated sequences of plastid protein-coding genes ( Fig. 7 and Supplementary Fig. S1). However, with regard to inter-generic relationships, some topological incongruences were found between cp genomes and nrDNA ITS sequences, which were in agreement with the earlier results on nrDNA ITS 25,26,28 and cpDNA marker 27,28 phylogenies. In both plastome based phylogenies, Reichardia was sister to Launaea, and the clade of Reichardia and Launaea was sister to Sonchus (Fig. 7a), while Reichardia was sister to the clade containing Launaea and Sonchus in nrDNA ITS phylogeny (Fig. 7b). It is likely difficult to decisively determine the inter-generic relationships in this study based on the current sampling, which included only two Launaea species. Given its species diversity (approximately 54 species and 10 subspecies), life form variability (perennial and annual herbs, subshrubs, cushion-forming rosette shrubs, and spinescent shrubs), and chromosome variation (n = 9, 8, 7, 6, and 5) 50 , the phylogenetic position of the genus Launaea remains to be clarified in future studies.
Species relationships within Reichardia also demonstrated incongruences between the cp genome and nrDNA ITS phylogenies. One novel finding of the current plastid phylogenomic study was that species with different chromosome numbers of n = 7 (R. picroides), n = 8 (R. intermedia), and n = 9 (R. albanica) were resolved in the same clade. This clade was sister to another clade with n = 8 species exclusively. In the nrDNA ITS phylogeny reconstructed previously 19 and in this study (Fig. 7b), a clade composed of species with n = 9 chromosomes (R. albanica, R. macrophylla, and R. dichotoma) was resolved as sister to a clade of the species with chromosome numbers of n = 8 and n = 7, clearly supporting the descending dysploidy process for chromosomal evolution in Reichardia. One of the derived groups included two species with n = 7 and 8 (R. picroides and R. intermedia, respectively), whereas the other was exclusively comprised of five n = 8 species (R. ligulata, R. famarae, R. crystallina, R. tingitana, and R. gaditana). However, the current study strongly suggested a close relationship among n = 7, n = 9, and a representative with n = 8 (R. intermedia) species on both plastome-based phylogenies ( Fig. 7a and Supplementary Fig. S1). According to our results, two major lineages may have diverged early from the common ancestor of the Reichardia species, that is, one lineage with n = 8 exclusively and the other lineage containing species with n = 9, 8, and 7. These relationships were not in agreement with a simple sequential descending dysploidy pattern from n = 9 to n = 7, as suggested by the nrDNA ITS phylogeny. The ancestral Reichardia karyotype was suggested to be equivocal, but n = 9 was most likely (56.5% probability at node 40; Fig. 8b and Supplementary Table S4) from the reconstruction mapped onto the nrDNA ITS ML tree. Even from the plastome tree, it was considered equivocal, but with the highest probability for n = 9 (27.4% at node P31), followed by n = 8 (24.6%), and n = 7 (17.3%). Based on these results, we can hypothesize that the common ancestor of Reichardia (with the most likely chromosome number of n = 9) initially diverged into two lineages: one with n = 8 species and the other including n = 9, n = 8, and n = 7 species. The hypothesis that n = 9 is the ancestral chromosome number in Reichardia coincides with the closely related genera Launaea and Sonchus, where this number is the most common and widely distributed 13,50,51 . Furthermore, it has been assumed as the ancestral basic chromosome number of the tribe Cichorieae 52-54 .
Given their current geographical distribution ranges, the affinity of n = 9 species with R. intermedia (n = 8) and R. picroides (n = 7) based on plastome sequences is an unexpected finding. The last two species show a wide distribution in large circum-Mediterranean regions in Africa, Asia, and Europe which hardly overlapped with those of n = 9 species and only partly with R. dichotoma in certain regions of West Asia (Turkey and Lebanon-Syria) [3][4][5] . In contrast to these two species, n = 9 species are currently found in very limited and isolated areas with a disjunctive geographical distribution pattern; Albanian endemic R. albanica is found only in Mount Çika in southern Albania at altitudes between 1000 and 1700 m above sea level (a.s.l.), Balkan endemic R. macrophylla in the central part of the Balkan Peninsula (East Dinarides) at altitudes between 500 and 1750 m a.s.l., and Western Asiatic R. dichotoma in Turkey, Lebanon-Syria, Iran, and North and Trans-Caucasus 7 . Considering the putative ancestral Reichardia karyotype (n = 9) and relic species restricted to the eastern Mediterranean 17 , it is highly plausible that Reichardia originated in the eastern Mediterranean. The mean divergence time based on the nrDNA ITS sequences for the two closely related genera Launaea and Sonchus was estimated to be 6.4 million years (myr) 55 . The estimated divergence time of Reichardia from its close relatives could be earlier than 6.4 myr based on the ITS tree topology (i.e., Reichardia is sister to the clade of Launaea and Sonchus) or younger than 6.4 myr on the plastome tree (i.e., Sonchus is sister to the clade of Launaea and Reichardia). Thus, according to the proposed divergence estimation, we can hypothesize that Reichardia originated around the Messinian salinity crisis in the Mediterranean basin 56 . It is plausible that after its origin in the eastern region, Reichardia rapidly expanded its range and diversified during dramatic climatic changes in the Mediterranean, forming n = 8 (R. intermedia) and n = 7 (R. picroides) lineages. In addition, the rapid speciation of n = 8 species in the western Mediterranean region, including species endemic to Macaronesian Islands (R. famarae, R. crystallina, and R. ligulata), most likely occurred after initial splitting from the n = 9 lineages. The Macaronesian Islands are well known for adaptive radiation and the diversification of explosive insular species of numerous plant lineages, including Aeonium 57 , Argyranthemum 58,59 , Echium 60,61 , Tolpis 62 and woody Sonchus alliance [25][26][27][28]  www.nature.com/scientificreports/ Reichardia relative to Launaea and Sonchus and n = 9 likely ancestral may also suggest that the common ancestor of Reichardia existed widely in the circum-Mediterranean region in the late Miocene and subsequently diversified during climatic turmoil, forming a relictual n = 9 lineage in the eastern Mediterranean and speciating in the western part of the Mediterranean. Furthermore, we cannot ignore the role of Quaternary climatic oscillations in shaping distribution areas, particularly in current relicts. The wide distribution range of R. tingitana, which roughly coincides with the paleo-geographical limits of the Mediterranean, seems to indicate, at least partially, that the expansion and diversification of this lineage may not have been recent and must have taken place prior to the present geologic configuration of the Mediterranean region. Given the phylogenetic position of R. tingitana (i.e., positioned in a recently derived clade), we cannot completely rule out the possibility of its recent rapid range expansion during climatic oscillation in the region.
The phylogenetic relationships among Reichardia species based on plastome sequences were also supported by morphological and chemotaxonomic data. Achene heteromorphy was less marked in the n = 9 species than in the other groups. R. intermedia and R. picroides showed plainly smooth inner achenes with differentiated shapes (Fig. 1c). In addition, species with purple base ligules were exclusively present in R. tingitana, R. gaditana, and the Canary Islands endemic clade. Recio et al. 24 divided Reichardia species into two groups by scanning electron microscopy, in support of the findings of the present study. In group I, which contained R. intermedia and R. picroides, the outer achene surface was made up of isodiametric papillated cells, while the cylindrical/ conic-truncated inner cells showed long and more or less rectangular cells without papillae. They also included R. macrophylla in group I, although the inner achenes did not show the same level of distinctiveness. Contrarily, group II contained n = 8 species; R. ligulata, R. famarae, R. crystallina, R. tingitana, and R. gaditana showed similar microscopic surface traits on outer and inner achenes, that is, more or less isodiametric cells, which were clearly distinguishable in shape, surface, and color ( Fig. 1c; fruits of R. tingitana). The separation of Reichardia species into two major groups seems to be corroborated further by differences in chemical compositions between R. picroides (group I; n = 7, 8, and 9 species) and R. tingitana (group II; n = 8 species) 24 .
Lastly, the phylogenetic incongruence between the plastome and nrDNA ITS phylogenies could be the result of past gene flow occurring in Reichardia species. Hybridization and introgression often appear to be responsible for the phylogenetic incongruence between organellar and nuclear markers, as observed in many groups of plants; however, other factors, such as phylogenetic sorting or the retention of ancestral polymorphisms, could also play a role in generating incongruences 63,64 . The R. intermedia and R. picroides group shared the most recent common ancestor with n = 9 species in plastome-based phylogenies with strong support values (BS 99% or 97% in ML and PP 1 in BI trees) ( Fig. 7a and Supplementary Fig. S1), whereas they were more closely related to the n = 8 group of species in the nrDNA ITS phylogeny (BS 98% in ML and PP 1 in BI trees) (Fig. 7b). Considering the wide distribution range of both groups in large circum-Mediterranean regions in Africa, Asia, and Europe, they may easily overlap in terms of distribution. Therefore, it is plausible that R. intermedia and R. picroides captured the chloroplast genome from the ancient n = 9 group with a much broader distribution, resulting in a phylogenetic incongruence between the cpDNA and ITS phylogeny. Several cases of hybridization among Reichardia species have been documented, further supporting the likelihood of this hypothesis, including R. × baetica Gallego & Talavera (R. tingitana × R. intermedia), R. x canariensis Gallego & Talavera (R. tingitana × R. ligulata), and R. × sventenia Gallego & Talavera (R. tingitana × R. famarae) 6 . While the R. × sventenia and R. × canariensis cases are examples of hybridization between widely distributed R. tingitana and the locally restricted species endemic to the Canary Islands, the hybrid origin of R. × baetica occurred between two divergent n = 8 lineages with a much wider geographical distribution.
In conclusion, we found highly conserved plastomes, including gene order and content, in Reichardia and Launaea. We achieved the first full resolution regarding inter and intra-generic phylogenetic relationships based on the complete plastome sequences of most Reichardia species with variable chromosome numbers (n = 7, 8, and 9). The plastid phylogenomics strongly suggested the early divergence of two major lineages of Reichardia, one with a group of species with exclusively n = 8 chromosome numbers and the other with n = 9, 8, and 7 species together. Although the plastome-based relationships were not in full agreement with the stepwise descending dysploid hypothesis in Reichardia, the ancestral Reichardia karyotype was highly likely to be n = 9 species. It is necessary to utilize whole nuclear genome data and a thorough synteny analysis to elucidate convincingly the species relationships and chromosomal evolution within the genus Reichardia. Based on thorough characterization and comparative analyses of the plastomes, we discovered several informative mutation hotspots, which will increase the efficiency and feasibility of phylogenetic reconstruction among Reichardia and closely related genera. www.nature.com/scientificreports/ The annotated GenBank (NCBI, Bethesda, MD, USA) format sequence file was used to draw a circular plastid genome map (Fig. 1a) using the OGDRAW software v1.2 (CHLOROBOX) 69 .

Comparative plastome analyses. We performed several comparative plastome analyses of the eight
Reichardia, representing variable chromosomal numbers. The analyses also included two related Launaea plastomes (L. arborescens and L. nudicaulis). The codon usage frequency was calculated using MEGA7 70 with the relative synonymous codon usage (RSCU) value, which is the relative frequency of occurrence of the synonymous codon for a specific amino acid. The online program predictive RNA editor for plants (PREP) suite 71 was used to predict the potential RNA editing sites for annotated protein-coding genes with 35 reference genes available with known edit sites, based on a cutoff value of 0.8 (suggested as optimal for PREP-Cp). The overall sequence divergence was estimated using the Lagan alignment mode 72 in mVISTA 39 . The nucleotide diversity (Pi) was calculated using sliding window analysis (window length = 1000 bp and step size = 200 bp, excluding sites with alignment gaps) to detect the most divergent regions (i.e., mutation hotspots) in DnaSP 38 . To evaluate the natural selection pressure in the protein-coding genes of the ten plastomes, site-specific models implemented in EasyCodeML 41 were used in the preset running mode based on CodeML algorithms 42 . Seven codon substitution models with heterogeneous ω values across sites were investigated and compared to detect positively selected sites based on likelihood ratio tests (M0, M1a, M2a, M3, M7, M8, and M8a).
Phylogenetic analyses and ancestral state reconstructions. The phylogenetic positions of the newly sequenced plastomes of Reichardia and Launaea assembled in this study were investigated in the context of their relationships with closely related species in the subtribe Hyoseridinae, including the outgroup species Lactuca sativa from the same tribe (Cichorieae). We further compared the plastome phylogeny with the nrDNA ITS phylogeny to explore phylogenetic incongruence and gain insights into chromosome evolution. To reconstruct the plastome phylogenies, we analyzed 19 plastomes of major genera of the subtribe Hyoseridinae and one outgroup based on whole plastome sequences, as well as concatenated sequences of 80 common protein-coding genes, including only nine newly sequenced species and downloaded plastomes from GenBank (Table 1 and  Supplementary Table S5). The nrDNA ITS phylogeny was also reconstructed using the retrieved sequences of 26 accessions from GenBank as specified in Supplementary Table S5, which were selected similarly as the taxa in the plastome phylogeny for the species belonging to the genera Reichardia, Launaea, and Sonchus s.l. in the subtribe Hyoseridinae. The plastome and nrDNA ITS sequences were aligned using MAFFT v. 7 73 , and phylogenetic trees were constructed using IQ-TREE v. The evolutionary model of GTR + I + Г (General Time Reversible substitution model with a proportion of invariable sites and a gamma-distributed rate across sites) was selected for Bayesian inference. The Markov chain Monte Carlo simulation (MCMC) had a length of 4,000,000 generations, sampling every 100 generations, for both plastome and nrDNA ITS phylogeny. The average standard deviation of the split frequencies reached below 0.01, which indicates that two runs reached stationarity. The first 25% of samples were discarded as burn-in, and the remaining samples were retained for the construction of a 50% majority-rule consensus tree with clade frequencies.
The reconstructed ML trees for the plastome and nrDNA ITS sequences were used as phylogenetic hypotheses to perform ancestral state reconstruction of chromosome numbers using the R toolkit 74 and MrBayes ancestral states with R (MBASR) 75 . The coded trait scores of chromosome numbers of each taxon were mapped onto both phylogenies, and the ancestral states and shared-derived karyotypes were inferred using MrBayes 3.2.7a 48 . The parameters in the analyses were set to "character.type = unordered" and "n.samples = 500".

Data availability
The plastome datasets sequenced and analyzed in the current study are available from GenBank under the accession numbers specified in Table 1.